# coding=utf-8
import arcpy
import os
import sys
import time


def get_system_encoding():
    """获取系统默认编码，优先使用UTF-8"""
    return 'utf-8'  # 强制使用UTF-8，避免GBK编码问题


# 设置系统编码
SYSTEM_ENCODING = get_system_encoding()
print("检测到系统编码: %s" % SYSTEM_ENCODING)

# 重新加载sys模块以设置默认编码
reload(sys)
sys.setdefaultencoding(SYSTEM_ENCODING)

# 设置控制台输出编码（Windows环境）
if sys.platform == 'win32':
    import ctypes

    try:
        # 设置控制台编码为UTF-8
        ctypes.windll.kernel32.SetConsoleOutputCP(65001)
        print("已设置控制台编码为UTF-8")
    except:
        print("控制台编码设置失败，使用UTF-8")


def safe_decode(s, encoding=None):
    """安全解码函数，自动处理编码问题"""
    if not encoding:
        encoding = SYSTEM_ENCODING
    if isinstance(s, unicode):
        return s
    try:
        return s.decode(encoding)
    except UnicodeDecodeError:
        try:
            return s.decode('utf-8')  # 尝试使用UTF-8解码
        except:
            return s.decode(encoding, 'ignore')  # 忽略无法解码的字符


def safe_encode(s, encoding=None):
    """安全编码函数，自动处理编码问题"""
    if not encoding:
        encoding = SYSTEM_ENCODING
    if isinstance(s, str):
        return s
    try:
        return s.encode(encoding)
    except UnicodeEncodeError:
        return s.encode(encoding, 'ignore')  # 忽略无法编码的字符


def get_feature_classes(database):
    """获取MDB数据库中的所有要素类名称"""
    feature_classes = []
    try:
        # 设置工作空间为MDB数据库
        arcpy.env.workspace = safe_encode(database)

        # 获取所有要素类
        for fc in arcpy.ListFeatureClasses():
            fc_name = safe_decode(fc)
            feature_classes.append(fc_name)

        print("在数据库 %s 中找到 %d 个要素类" % (safe_decode(database), len(feature_classes)))
        return feature_classes
    except Exception as e:
        print("获取要素类列表时出错: %s" % safe_decode(str(e)))
        return []


def compare_same_name_fc(db_old, db_new, output_db, tolerance):
    """对比两个MDB数据库中相同名称的要素类，提取变化要素"""
    try:
        # 确保输出数据库存在
        if not arcpy.Exists(safe_encode(output_db)):
            arcpy.CreateFileGDB_management(os.path.dirname(safe_encode(output_db)),
                                           os.path.basename(safe_encode(output_db)))
        print("开始对比分析...")
        start_time = time.time()

        # 获取两个数据库中的要素类列表
        fc_list_old = get_feature_classes(db_old)
        fc_list_new = get_feature_classes(db_new)

        # 找出同名要素类
        common_fcs = set(fc_list_old).intersection(set(fc_list_new))

        if not common_fcs:
            raise ValueError("两个数据库中未找到同名要素类")

        total_fcs = len(common_fcs)
        print("找到 %d 个同名要素类，开始处理..." % total_fcs)

        # 遍历处理每个同名要素类
        completed_count = 0
        for fc_name in common_fcs:
            process_feature_class(db_old, db_new, fc_name, output_db, tolerance)
            completed_count += 1
            print("已处理: %d/%d 个要素类" % (completed_count, total_fcs))

        elapsed_time = time.time() - start_time
        print("\n所有要素类对比完成，共处理 %d 个要素类，总耗时: %.2f 秒" % (total_fcs, elapsed_time))

    except arcpy.ExecuteError:
        error_msg = arcpy.GetMessages(2)
        print_error(error_msg)
    except Exception as e:
        error_msg = str(e)
        print_error(error_msg)


def print_error(msg):
    """安全打印错误信息，处理编码问题"""
    try:
        # 先尝试用UTF-8解码
        safe_msg = safe_decode(msg, 'utf-8')
        print("处理过程中发生错误: %s" % safe_msg)
    except:
        # 再尝试用系统编码解码并忽略错误
        safe_msg = safe_decode(msg, 'ignore')
        print("处理过程中发生错误: %s" % safe_msg)


def process_feature_class(db_old, db_new, fc_name, output_db, tolerance):
    """处理单个同名要素类的对比，提取并合并变化要素"""
    try:
        # 确保要素类名不含非法字符
        clean_fc_name = ''.join(c if c.isalnum() or c == '_' else '_' for c in fc_name)
        fc_name_encoded = safe_encode(clean_fc_name)
        print("\n===== 处理要素类: %s =====" % fc_name_encoded)

        # 构建要素类路径
        fc_old = os.path.join(safe_encode(db_old), fc_name_encoded)
        fc_new = os.path.join(safe_encode(db_new), fc_name_encoded)
        output_change = os.path.join(safe_encode(output_db), "%s_change" % fc_name_encoded)  # 变化要素输出路径

        # 检查要素类是否存在
        if not arcpy.Exists(fc_old) or not arcpy.Exists(fc_new):
            print("警告: %s 在任一数据库中不存在，跳过处理" % fc_name_encoded)
            return

        # 获取几何类型并正确显示中文
        desc_old = arcpy.Describe(fc_old)
        geom_type = desc_old.shapeType
        geom_type_cn = {
            "Point": "点",
            "Multipoint": "多点",
            "Polyline": "线",
            "Polygon": "面"
        }.get(geom_type, geom_type)
        print("%s 几何类型: %s" % (fc_name_encoded, geom_type_cn))

        # 检查空间参考是否一致
        desc_new = arcpy.Describe(fc_new)
        if desc_old.spatialReference.factoryCode != desc_new.spatialReference.factoryCode:
            print("警告: %s 空间参考不一致，可能影响对比结果" % fc_name_encoded)

        # 获取字段列表
        fields = [safe_decode(f.name) for f in arcpy.ListFields(fc_new) if not f.required]
        fields.insert(0, "SHAPE@")
        fields = [safe_encode(f) for f in fields]

        # 提取新增和消失要素
        added_features = extract_added_features(fc_old, fc_new, fields, tolerance, geom_type)
        deleted_features = extract_deleted_features(fc_old, fc_new, fields, tolerance, geom_type)

        # 合并新增和消失要素
        all_change_features = added_features + deleted_features
        change_count = len(all_change_features)

        if change_count == 0:
            print("警告: 未检测到变化要素，跳过创建文件")
            return

        # 创建输出要素类（变化要素）
        arcpy.management.CreateFeatureclass(
            safe_encode(output_db),
            os.path.basename(safe_encode(output_change)),
            geometry_type=geom_type,
            template=fc_new,
            spatial_reference=desc_new.spatialReference
        )

        # 保存变化要素
        merge_features(all_change_features, output_change, fields, geom_type)

    except Exception as e:
        error_msg = str(e)
        print("处理 %s 时出错: %s" % (fc_name_encoded, safe_decode(error_msg, 'utf-8')))


def extract_added_features(fc_old, fc_new, fields, tolerance, geom_type):
    """提取新增要素"""
    added_features = []

    if geom_type == "Point":
        # 创建旧数据的点坐标集合（带容差）
        point_coords_old = set()
        with arcpy.da.SearchCursor(fc_old, ["SHAPE@XY"]) as cursor:
            for row in cursor:
                x, y = row[0]
                x_tol = round(x / tolerance) * tolerance
                y_tol = round(y / tolerance) * tolerance
                point_coords_old.add((x_tol, y_tol))

        # 提取新增点
        with arcpy.da.SearchCursor(fc_new, fields) as search_cursor:
            for row in search_cursor:
                geom = row[0]
                x, y = geom.firstPoint.X, geom.firstPoint.Y
                x_tol = round(x / tolerance) * tolerance
                y_tol = round(y / tolerance) * tolerance

                if (x_tol, y_tol) not in point_coords_old:
                    added_features.append(row)

    else:
        # 合并旧数据的几何为单个要素
        union_geom_old = merge_geometries(fc_old, geom_type)
        if not union_geom_old:
            union_geom_old = arcpy.Geometry()  # 创建空几何

        # 提取新增要素
        with arcpy.da.SearchCursor(fc_new, fields) as search_cursor:
            for row in search_cursor:
                geom = row[0]
                if not geom.overlaps(union_geom_old) and not geom.touches(union_geom_old):
                    added_features.append(row)

    print("提取到 %d 个新增要素" % len(added_features))
    return added_features


def extract_deleted_features(fc_old, fc_new, fields, tolerance, geom_type):
    """提取消失要素"""
    deleted_features = []

    if geom_type == "Point":
        # 创建新数据的点坐标集合（带容差）
        point_coords_new = set()
        with arcpy.da.SearchCursor(fc_new, ["SHAPE@XY"]) as cursor:
            for row in cursor:
                x, y = row[0]
                x_tol = round(x / tolerance) * tolerance
                y_tol = round(y / tolerance) * tolerance
                point_coords_new.add((x_tol, y_tol))

        # 提取消失点
        with arcpy.da.SearchCursor(fc_old, fields) as search_cursor:
            for row in search_cursor:
                geom = row[0]
                x, y = geom.firstPoint.X, geom.firstPoint.Y
                x_tol = round(x / tolerance) * tolerance
                y_tol = round(y / tolerance) * tolerance

                if (x_tol, y_tol) not in point_coords_new:
                    deleted_features.append(row)

    else:
        # 合并新数据的几何为单个要素
        union_geom_new = merge_geometries(fc_new, geom_type)
        if not union_geom_new:
            union_geom_new = arcpy.Geometry()  # 创建空几何

        # 提取消失要素
        with arcpy.da.SearchCursor(fc_old, fields) as search_cursor:
            for row in search_cursor:
                geom = row[0]
                if not geom.overlaps(union_geom_new) and not geom.touches(union_geom_new):
                    deleted_features.append(row)

    print("提取到 %d 个消失要素" % len(deleted_features))
    return deleted_features


def merge_features(change_features, output_fc, fields, geom_type):
    """合并变化要素到输出要素类"""
    change_count = len(change_features)
    print("检测到 %d 个变化要素，开始保存..." % change_count)

    # 写入合并后的要素
    with arcpy.da.InsertCursor(output_fc, fields) as insert_cursor:
        for row in change_features:
            insert_cursor.insertRow(row)

    print("已将 %d 个变化要素保存至: %s" % (change_count, safe_encode(output_fc)))


def merge_geometries(fc, geom_type):
    """合并要素类的所有几何为单个几何"""
    try:
        geometry_list = []

        # 读取所有几何
        with arcpy.da.SearchCursor(fc, ["SHAPE@"]) as cursor:
            for row in cursor:
                geometry_list.append(row[0])

        if not geometry_list:
            print("警告: 数据为空，无法合并")
            return None

        # 创建临时要素类
        temp_fc = os.path.join(arcpy.env.scratchGDB, "temp_merge")
        arcpy.CreateFeatureclass_management(
            arcpy.env.scratchGDB,
            "temp_merge",
            geom_type,
            spatial_reference=geometry_list[0].spatialReference
        )

        # 写入几何
        with arcpy.da.InsertCursor(temp_fc, ["SHAPE@"]) as icursor:
            for geom in geometry_list:
                icursor.insertRow([geom])

        # 根据几何类型选择合并工具
        union_output = os.path.join(arcpy.env.scratchGDB, "union_result")
        if geom_type in ["Polygon", "Polyline"]:
            arcpy.Union_analysis(temp_fc, union_output, "ALL")
        elif geom_type == "Point":
            arcpy.Dissolve_management(temp_fc, union_output, None, None, "SINGLE_PART", "DISSOLVE_LINES")
        else:  # Multipoint
            arcpy.CopyFeatures_management(temp_fc, union_output)

        # 读取合并后的几何
        with arcpy.da.SearchCursor(union_output, ["SHAPE@"]) as cursor:
            union_geom = cursor.next()[0]

        # 清理临时数据
        arcpy.Delete_management(temp_fc)
        arcpy.Delete_management(union_output)

        return union_geom

    except Exception as e:
        error_msg = str(e)
        print("几何合并失败: %s" % safe_decode(error_msg, 'utf-8'))
        return None


if __name__ == "__main__":
    # 配置参数（MDB数据库路径）
    db_19 = r"D:\ArcGis\DATA\19_data.mdb"  # 19年数据MDB数据库
    db_24 = r"D:\ArcGis\DATA\24_data.mdb"  # 24年数据MDB数据库
    output_db = r"D:\ArcGis\Output\change_result.mdb"  # 输出结果MDB数据库
    TOLERANCE = 0.000000001  # 点要素对比容差（米）

    # 执行对比分析
    compare_same_name_fc(db_19, db_24, output_db, TOLERANCE)